Generalized Dynamical Mean-Field Theory for Bose-Fermi Mixtures 

in Optical Lattices 
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We give a detailed discussion of the recently developed Generalized Dynamical Mean-Field Theory 
(GDMFT) for a mixture of bosonic and fermionic particles. We show that this method is non- 
perturbative and exact in infinite dimensions and reliably describes the full range from weak to 
strong coupling. Like in conventional Dynamical Mean-Field Theory, the small parameter is 1/z, 
where z is the lattice coordination number. We apply the GDMFT scheme to a mixture of spinless 
fermions and bosons in an optical lattice. We investigate the possibility of a supersolid phase, 
focusing on the case of 1/2 filling for the fermions and 3/2 filling for the bosons. 
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I. INTRODUCTION 

The impressive experimental progress in the field of 
ultracold atoms in the last decade has brought it to 
the forefront of research on strongly correlated quan- 
tum many-body systems. The possibility to confine and 
manipulate atoms in optical lattices created by stand- 
ing waves of laser light gives the opportunity to real- 
ize some of the model Hamiltonians of condensed matter 
physics, and in this way shed light on notoriously difficult 
problems^ 2 -" 3 -. Going beyond that, also systems without 
clear analogue in condensed matter systems can be real- 
ized. 

A prime example of this is the possibility to study 
bosonic atoms in an optical lattice 1 -' 2 -' 4 -'^' 7 -'^. These 
systems allow for the experimental check of theoretical 
predictions on the Bose-Hubbard model, introduced by 
Fisher et al3 in the late eighties. Recent experiments 
with high accuracy verified the supcrfluid-Mott insula- 
tor phase transition 2 ^ 5 -. These experimental results are 
in good agreement with theoretical predictions for the 
Bose-Hubbard mode&i&ii. 

Cold atomic gases also offer the possi- 
bility to realize mixtures of fermions and 
bocon a 12 - 13,14 - 15,16 - 17,18 - 19,20 - 21,22 - 23,24 - 25 - 26,27 This 
yields a very rich system, which at this moment is 
far from fully explored. One promising route that 
is currently experimentally investigated is to make 
heteronuclear molecules consisting of a boson and a 
fermion, with a permanent electrical dipole moment and 
hence a long range interaction^ 2 .. In this paper we will, 
however, concentrate on the many-body behavior of an 
interacting cloud of spinless fermions and bosons. 

This system bears some analogy with the well- 
known two-component Fermi-Fermi mixture, but is in 
fact much richer. By replacing one of the fermionic 
components by bosons, one keeps the instability of 
half-filled fermions towards charge-density wave (CDW) 
ordering. For historical reasons we keep this termi- 
nology throughout this paper, although the fermionic 
atoms under consideration do not carry a charge. At 
the same time the bosonic species can be supcrfiuid, 



allowing for supersolid behavior, where diagonal CDW 
order coexists with off-diagonal superfluid long-range 
order. Several previous theoretical works have stud- 
ied mixtures of fermions and bosons in an optical 

1 a ^^ c o 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 4t. 42. 43.44. 45. 46.47 

In some of thes o 34 ' 35 i 36 ' 37 i 38 supersolid phases were 
predicted. 

Investigating a strongly correlated Bosc- Fermi mixture 
in an optical lattice is a difficult problem, to which pow- 
erful numerical and analytical techniques have been ap- 
plied. In one dimension this involved bosonization 3 ^, 
Density Matrix Renormalization Group 30 i 33 , and Quan- 
tum Monte Carl o 38 ' 39 ! 40 ' 41 ^ 42 ' 43 . In higher dimensions, 
however, non-perturbative calculations are sparse. In 
two dimensions Renormalization Group studie a 46 i 47 have 
been carried out. Although able to describe non- 
perturbative effects, this technique is limited to weak 
couplings. Another powerful technique that has been ap- 
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plied in twos 7 -, and recently also three dimension; 
to integrate out the fermions. In this way one generates 
a long-ranged, retarded interaction between the bosons, 
which means that the resulting bosonic problem is still 
hard to solve. Important progress has recently been made 
in mapping out the Mott-insulating lobes. A compos- 
ite fermion approach^ 5 , was used to qualitatively describe 
possible quantum phases of the Bose-Fermi mixture. 

In this paper we describe the recently introduced Gen- 
eralized Dynamical Mean-Field Theory (GDMFT^i to 
study this system. This is a non-perturbative method 
which becomes exact in infinite dimensions and is a good 
approximation for three spatial dimensions. The only 
small parameter is 1/z, where z is the coordination num- 
ber. For this reason, the method reliably describes the 
full range from weak to strong coupling. To solve the ef- 
fective self-consistent quantum impurity problem arising 
within GDMFT, we use the Numerical Renormalization 
Group (NRG)^ 8 -. NRG resolves the low-frequency infor- 
mation very well, which enables us to reliably capture 
the supersolid phase, which in general has a small gap. 

The paper is organized as follows: in the next section 
we will shortly describe the Hamiltonian of the system 
and afterwards in section ILnl we consider GDMFT in de- 



tail. In Sec. [TVl wc apply the GDMFT to a mixture of 
spinlcss fermions and bosons at commensurate filling, in 
particular for the case when the fermions arc half-filled, 
while the filling of the bosons is 3/2. In Sec. [V]we end 
up with concluding remarks. In Appendix fA] we derive 
the effective action, while in Appendix [B] and Appendix 
[Cl we derive the expression for the kinetic energy and 
self-energy, respectively. 



II. MICROSCOPIC MODEL 

The standing waves of an optical lattice produce a po- 
tential V b ^(r) = V^ f \sin 2 (kx) + sin 2 (ky) + sin 2 (kz)), 
with k = 27r/A where A is the wavelength of the laser. 
Throughout this paper we assume the optical lattice to be 
strong enough that we can restrict ourselves to the low- 
est band. This means that we require Vq' /E^' > 2, 

where E^' = H 2 k 2 jlm^f) is the recoil energy for bosons 
(fermions). In order for the single band approximation to 
hold, all the other energy scales and temperatures should 
be smaller than the band gap. Since the Wannier func- 
tions for the fermions and the bosons arc well localized, 
it is a good approximation to consider only local interac- 
tions between particles and next-neighbor hopping, i.e., 
to consider the system in a tight-binding approximation. 
Under these approximations a mixture of fermions and 
bosons in an optical lattice can be well described by the 
single-band Fcrmi-Bosc Hubbard model 

W = - E { t f4^ + tbWb j }-J2{^ f n{ + ^} 

+ E{f^ ( ^ _ 1} + u f*ti h ii + r } ' (1) 

where denotes summation over nearest neighbors. 

cj CT (fe|) is the fermionic (bosonic) creation operator at 
site i, while n{ a = c\ a Ci a {n\ = b\bi) denotes the number 
operator for fermions and bosons, and n{ = + is 
the total fermionic particle number on site i. fj,b and 
denote the chemical potentials for bosons and fermions, 
respectively. Ub, Uf and E//& are the on-site boson-boson, 
fcrmion-fcrmion and fcrmion-boson interactions, respec- 
tively and tf(h\ is the tunneling amplitude for fermions 
(bosons). The following relation holds between the pa- 
rameters of the model and the experimental parameters: 
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FIG. 1: (Color online) Schematic picture of Generalized Dy- 
namical Mean-Field theory (GDMFT): within the GDMFT 
approach the full many-body lattice problem is replaced by a 
single-site problem, which is coupled to the fermionic bath as 
in "usual" DMFT and to the bosonic bath via the Gutzwiller 
approach. 

where a^, a/ and a/& are boson-boson, fermion-fermion 
and fermion-boson scattering lengths. From Eqs. it 
is clear that the ratio of the interaction to the tunneling 
amplitude can be varied from weak to strong coupling. 

In the case of spinless fermions, since there is only one 
species of fermions and the interaction is purely local, 
the fermionic part simply reduces to the free fermionic 
Hamiltonian. The total Hamiltonian therefore has the 
following form 

+E{y^ 6 - 1 )+ c/ /^}' ( 5 ) 



III. METHOD 

A. Self-consistency loop 

Following the very successful Dynamical Mean-Field 
Theory (DMFT)i^ and Gutzwilleri£ schemes, which 
are exact in infinite dimensions, we consider first the 
infinite-dimensional limit {d — > oo) of the Bose-Fermi 
mixture, which is expected to be a good approximation 
to three spatial dimensions. The main idea of the DMFT 
approach is to map the quantum lattice problem with 
many degrees of freedom onto a single site - the "impu- 
rity site" - coupled self-consistently to a non- interacting 
bath. To derive the self-consistency equations for this 
model, we use the "cavity method" 49 ' 50 : one considers a 
single site of the lattice and integrates out the remaining 
degrees of freedom on all other sites. To derive the self- 
consistency relations, we use the path integral formalism. 
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The important point in this derivation is that we consider 
the limit of infinite spatial dimensions (i.e. lattice coordi- 
nation number z — > oo). To keep the kinetic energy finite, 
we need to rescale the hopping parameters of the Hamil- 
tonian JI) as follows: t f = t)/jz 4 a and t b = t* b /z 
Doing so, the parameter l/z appears as a small parame- 
ter in the theory, which is used to control the expansion. 
We note here that l/z is not a coupling parameter in the 
original Hamiltonian. Therefore this method is suited for 
the full range of couplings considered. This gives us also 
a way to estimate accuracy: neglecting terms of order l/z 
leads to reasonably small errors for the three-dimensional 
cubic lattice where z = 6. 

The first step in this formalism is to derive the effective 
action of the impurity site (for details see Appendix lX]) by 
integrating out the remaining degrees of freedom (i ^ 0) 
in the partition function: 



Seff - 



Jeff 



1/ II DclDc^DbtDbi 



e- s . (6) 



i=£0,o 



where Cj CT , c* CT are Grassmann variables describing 
fermions, bi, b* are C-numbers describing bosons. To 
leading order in l/z one obtains 

S eff = _ H/ dri dT2 Yl ^oA^G^in - T 2 )c 0a (T 2 ) 
a J JO 

+ / dr bo{T)(d T - Hb)b (T) 
Jo 

-t b [ drY!^i(rMr)+c.c)+Uf[ drn^ (r)n^ (r) 
Jo ■ Jo 

+U fb [ dr4{r)n\{r) + U b /"drng(r)(ng(r) - 1). (7) 
Jo Jo 

Here 5>°(t) = (6)° is the bosonic superfluid parameter, 
which is static. We have introduced the Weiss function 

where G° iji<7 (n - r 2 ) = -{Tc ia {T X )c] a {T 2 ))° is the inter- 

acting Green's function for the fermions, and means 
summation only over the nearest neighbors of the "im- 
purity site" . The expectation values are here calculated 
in the cavity system without the impurity site, which is 
indicated by the notation (. . .)°. 

The next step in the derivation is that the expecta- 
tion values in the cavity system are identified with the 
expectation values on the impurity site. This means that 
we identify $?(r) = (b)° = (b) and G° ii>a (n - r 2 ) = 

-(Tc^rOCfc)} = -(Tc 0CT (r 1 )4(r 2 )) , where the 
notation (. . .)o means expectation value for the impurity 
site. In passing by, we note that this involves again an 
error of order l/z (vanishing in the limit of high dimen- 
sionality), since a site at the edge of the cavity has one 
neighbor less compared to the impurity site. However, 
in this way, we have derived a self-consistency relation, 
which only involves the impurity site. 




FIG. 2: (Color online) Schematic structure of the Bethe lat- 
tice (here with coordination number 2 = 3). 



By inspecting these self-consistency relations, it be- 
comes clear that the bosonic part corresponds to the 
Gutzwiller approximation, whereas the fermionic part 
corresponds to DMFT. The two are coupled by the on- 
site density- density interaction. We note here that this 
derivation shows that the Gutzwiller approximation for 
bosons is exact in infinite dimensions, and, like DMFT, 
valid for arbitrary couplings in the Hamiltonian. There- 
fore this approximation is able to describe the whole 
phase-diagram, in particular the transition from super- 
fluid to Mott-insulator. This point is not always appreci- 
ated in the literature, where the Gutzwiller approxima- 
tion is sometimes regarded as a strong-coupling approx- 
imation. 

Summarizing, the GDMFT employed in our calcula- 
tion consists of the DMFT algorithm for the fermions, 
combined with bosonic Gutzwiller mean-field theory. 
The bosons are described by the superfluid order param- 
eter 3>°(t) = (&(t)) while the fermions are characterized 
by the Weiss Green's function 

G~ l {iu; n ) = iuj n + \x a - t^y^ G° jt0 .(iu) n ) . (8) 



where oj n = (2n + 1)tt//3 arc Matsubara frequencies. 
Schematically the GDMFT is depicted in Fig. [TJ 

The self-consistency equation for the fermions assumes 
the simplest form for the Bethe lattice which is schemat- 
ically depicted in Fig. [5] and has a semi-elliptic non- 
interacting density of states p(e) = \J ^tj 2 — e 2 /2irtj 2 . 
The reason for this simplification is that for the Bethe 
lattice the summation in Eq. © is reduced to i = j, 
because all neighbors of "impurity site" are decoupled. 
The self-consistency relation for fermions on the Bethe 
lattice is therefore 



Q^iiWn) = iu}„ + [la - t* f G a (iU} n ) . 



(9) 



The self-consistent GDMFT loop has the following 
structure: we start from an initial guess of the Weiss 
Green's function and superfluid order parameter. The 
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effective action of the model is then given by Eq. (J7J), 
which allows us to calculate all local Green's functions 
and expectation values, including the interacting Green's 
function and the superfluid order parameter. The loop 
is closed by Eq. ([9]), from which we calculate the new 
Weiss Green's function. This procedure is repeated until 
convergence is reached. 

B. Generalized single impurity Anderson Model 

The most difficult step in the procedure outlined above 
is the calculation of the local Green's function from the 
effective action. We use the Numerical Renormalization 
Group for this purpose, which is non-perturbative and 
provides reliable low-frequency information. 

To be able to employ NRG, we map the self-consistent 
single-site model onto a Generalized Single Impurity An- 
derson Model (GSIAM), which by construction has ex- 
actly the same effective action ([7]) as the initial Hamil- 
tonian {!]). As in the conventional Single Impurity An- 
derson model (SIAM), the impurity site is coupled to a 
non-interacting fcrmionic bath which - like the effective 
action {?]) - needs to be determined self-consistently in 
Dynamical Mean-Field Theory. In addition, the GSIAM 
now also contains a bosonic degree of freedom on the "im- 
purity site" , which is self-consistently coupled to the su- 
perfluid order parameter, according to Gutzwillcr mean- 
field theory. In summary, the GSIAM is described by the 
following Hamiltonian, which allows for a two-sublattice 
structure: 

Wgsiam = ^ [Hb + Wfb + 

a=±l 

HZ = -zt b dp a bi + <p* s b a ) + - !) - w£ 

U% = U fb n{h b a 

n a f = -^ f h{ + U f n\ a h{ a + 

+ ^2{eia a a\ aa a llja + Vi aa (cl a a laa + h.c)j } 

Here a = ±1 is the sublattice index (a = —a), z is the 
lattice coordination number, ip a = (b a ) is the superfluid 
order parameter on sublattice a. I labels the noninter- 
acting orbitals of the effective bath, and Vi aa are the cor- 
responding fcrmionic hybridization matrix elements^. 

C. Numerical Renormalization Group 

The Hamiltonian Eq. (fT0|) can be diagonalized using 
the Numerical Renormalization Group (NRG)^. The 
key idea of this method is to perform a logarithmic dis- 
cretization of the conduction band in order to exploit the 
separation of energy scales crucial for a renormalization 
group treatment. By an additional unitary transforma- 
tion the conduction band is mapped onto a semi-infinite 



linear chain. The fcrmionic part of our Generalized An- 
derson Impurity Model (GSIAM) in Eq. (TO]) then takes 
the form 

UJ = -^ f n{ + U f n{ a h{ a + (11) 

+ e rLaa(yd\ l _ laa d niJa + h.cj + ^naa^naadnaa 

where d\ laa and e naa are the fermion creation operators 
and the hopping coefficients on the linear chain. ($_ laa = 
c\ a corresponds to the impurity site. S ntTa is the on- 
site energy for site n of the linear chain. Due to the 
logarithmic discretization, the hopping parameters and 
onsite energies now decay exponentially e naa ~ A~"/ 2 
and S naa ~ A~™/ 2 , where A is the NRG discretization 
parameter, which in our calculations we have chosen as 
A = 2. 

As is obvious from Eq. 1101 the bosons are incorporated 
only on the "impurity site" and self-consistently coupled 
to the superfluid order parameter. This means that they 
will not affect the renormalization scheme of NRG, but 
only the construction of the "impurity" Hamiltonian. In 
order to keep the dimension of the impurity site Hilbcrt 
space small enough to handle it numerically, we use a 
cut-off for the number of bosons on the impurity site, 
which can be kept low due to the repulsive interactions, 
which suppress multiple occupancy of the bosons. 

The renormalization scheme of NRG then works as 
follows^: In each step one more site of the linear chain 
is added to the Hamiltonian, and using the eigenvalues 
and the eigenvectors of the Hamiltonian in the previous 
step one can build the Hamiltonian for this system. The 
next step is to diagonalize the new Hamiltonian and find 
its eigenvalues and eigenvectors. The size of the Hilbert 
space after adding one more site increases by a factor 4 
for two-component fermions and by a factor 2 for spinlcss 
ones. To limit the matrix size, one then keeps only the 
Ni eve i lowest energy levels (usually Ni eve i = 600 — 1000) 
in each step. This truncation scheme is controlled by the 
energy scale separation discussed above. The number of 
iterations Ni ter is directly related to the temperature of 
the system as k B T ~ DA~ Nu<!r / 2 where D = 2t* denotes 
the fermionic half-bandwidth. In zero temperature calcu- 
lations, such as in this work, Ni ter is chosen large enough 
to yield a temperature below any intrinsic energy scale 
of the system. Here we have chosen Nu er = 60. 

From the eigenstates and matrix elements thus ob- 
tained one can then calculate any local expectation value 
or correlation function, such as the superfluid order pa- 
rameter tp a = (b a ) and the local fermionic interacting 
(impurity) spectral function A aa (uj) . 

D. Ground state energy 

It is clear that the final result of the GDMFT calcula- 
tions should not depend on the initial conditions of the 
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self-consistency loop. However, for physical reasons it 
can happen that the self-consistent GDMFT procedure 
yields multiple stable solutions. To find the ground state 
of the system in those cases, we need to compare the en- 
ergies of the coexisting solutions. The ground state will 
correspond to the solution with the lowest energy. For 
this purpose we need to calculate the total energy which 
is given as follows: 



N = ~7V 



(12) 



= -ztb(f-i(px + Ids ep(e) du)B a (e,uj) 

_ + \J — oc J — oc 



a=±l 

where the index a = ±1 corresponds to the two different 
sublattices. To calculate the fermionic part of the kinetic 
energy above, we have used the same approach as for an 
antifcrromagnetic state, which also has a two-sublattice 
struct ur o 50 i 53 (for details see Appendix [Bj . p(e) is the 
fermionic non-interacting density of states and 



B a {e, to) = Im 



1 



y/ Cr,lCo\-l - £ 



(13) 



is a spectral function, with ( aa = lu + fi a f — T, aa (uj). 
We calculate the self-energy as follows^(for details see 
Appendix [Cj : 



G aa {u) 



(14) 



where G aa {uj) = (f aa f% a )uj is a local fermionic single- 
particle Green's function, F/^(w) = (f aa fLfsafD^ 
and Fl b a (u) = (UJtbJlJw Here a = -a denotes the 
opposite spin state. 

For nonzero temperature (not considered here) the free 
energy is the relevant quantity, which means that also the 
entropy has to be calculated. 

E. Evaluation 

We close this section with a short summary of the 
method. The GDMFT technique is a combination of 
the DMFT and Gutzwiller approaches. We have shown 
that it is exact in infinite dimensions, and it is assumed 
to be a good approximation for three spatial dimen- 
sions (with the lattice coordination number z = 6). 
Fermionic DMFT calculations in three dimensions show 
indeed excellent agreement with QMC calculations^ and 
experiments^. The only small parameter in this method 
is 1/z (where z is the lattice coordination number). 
GDMFT therefore incorporates local correlations be- 
tween bosons and fcrmions in a fully non-perturbativc 
fashion. Non-local correlations, on the other hand, can 
be calculated only on a mean-field level. 



Since the fermions are treated with a dynamical mean- 
field, their quantum fluctuations are also captured. 
Higher orders in 1/z could make quantitative changes, 
but no qualitative changes are expected. The bosons on 
the other hand are treated in static mean field and cou- 
ple only to the bosonic order parameter. Although this 
is indeed exact in infinite dimensions, for a finite number 
of spatial dimensions even normal (i.e. non-supcrfluid) 
bosons will hop. This will e.g. affect the fluctuations in 
the boson number (ri?) — (h^) 2 . Within the Gutzwiller 
approximation this quantity is zero in the Mott insula- 
tor and alternating Mott insulator phase (which will be 
defined in the next section). The inclusion of normal 
hopping would lead to finite fluctuations. This effect is 
however not essential for the physics of the supersolid dis- 
cussed here. In future calculations, normal bosonic hop- 
ping could be included via the recently developed Bosonic 
DMFT (BDMFT)£I£I. 

The above derivation was valid independently of tem- 
perature and impurity solver. Therefore, GDMFT also 
gives a reliable description of a Bosc-Fcrmi mixture in an 
optical lattice at any finite temperature. As an impurity 
solver one can use NRG^ or exact diagonalizatio n 50 i 59 i 60 
which works very convenient at finite temperature. In the 
present work, we only apply it at T = and using NRG 
as an impurity solver. 



IV. SUPERSOLID AND ALTERNATING MOTT 
INSULATOR FOR 3/2-FILLED BOSONS 

A. GDMFT analysis 

The supersolid phase - the phase with coexisting bro- 
ken U(l) symmetry and particle wave density order - 
is one of the intriguing subjects in condensed matter 
physics. It is still an open question whether a super- 
solid has been realized in recent experiments on He 
61 . While in single-component quantum gases super- 
solids can only be stabilized by including nearest neigh- 
bor repulsion between the particles^, they can be con- 
veniently realized in Bose-Fermi mixtures with on-site 
repulsion in an optical lattice where the fermions are 
at half n i]i ng 34 1 35 1 36 1 3Z,38 1 4^ > The Hamiltonian for this 

mixture of bosons and spinless fermions is given in Eq. [5] 
The mechanism for supersolid formation here is the in- 
stability of fermions at half-filling towards charge-density 
wave (CDW) formation because of Fermi surface nest- 
ing. The bosons act as impurities for the fermions, which 
drives the system into this phase with broken transla- 
tional symmetry. Since the bosons remain supcrfluid 
for moderate interactions, the associated U{1) symme- 
try and the translational symmetry are simultaneously 
broken. For strong Bose-Fermi interactions, on the other 
hand, fermions and bosons avoid each other and are local- 
ized in different sublattices, thus forming an Alternating 
Mott Insulator (AMI) phase as shown before^. 

In our previous work^ we studied the Bose-Fermi 




FIG. 3: (Color online) Dependence of the amplitude of the 
bosonic/fermionic density wave on the Fermi-Bose interac- 
tion, for the case when ztb = 0.05D and Ub = 1.0D, where 
D denotes the half-band width of the fermions. In the inset 
we depict the superfluid order parameter. The different line 
types in the inset correspond to results on the two sublattices. 
The different phases are schematically depicted in Fig. [4] In 
this and all following plots energies are expressed in units of 
D. 

mixture for the case when both species were half-filled. 
We obtained three different phases: (i) Supersolid phase 
for small Bose-Fcrmi interaction and strong boson-boson 
interaction, (ii) AMI phase for strong Fermi-Bose and 
boson-boson interaction and (iii) phase separation for 
small boson-boson interaction. 

We remark here that those results, and also the results 
obtained in this paper, are obtained with a density of 
states without Van Hove singularities. In fact, the results 
were obtained using the density of states of the Bethe lat- 
tice, which is semi-elliptic and regular everywhere. We 
were still able to identify a supersolid phase, proving the 
point that a singularity in the non-interacting states is 
not a necessary condition for the occurrence of a super- 
solid. However, because of the lack of singularities in the 
density of states, the particle density oscillation and the 
gap in the spectrum in the supersolid phase were rather 
small. 

Therefore, in this paper we study a different case where 
the filling of fermions is 1/2, while the filling of the bosons 
is higher, namely (n£) = 3/2. The reason for this par- 
ticular choice is that it allows for two different Alter- 
nating Mott Insulator (AMI) phases, with amplitude of 
the bosonic density oscillation 1/2 and 3/2, respectively. 
These two AMI phases are separated by a supersolid 
phase. The amplitude of the density oscillations in this 
supersolid phase in between the two AMI phases is of 
order one, which makes the experimental detection much 
easier. 

We study the system using GDMFT— . To overcome 
the tendency towards phase separation in the system, 
we consider the case where the bosons are much slower 
than the fermions zt b = 0.05-D, and where the repulsion 
among the bosons is strong Ub = D. In Fig. Owe plot the 
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Supersolid-1 AMI-1 

WWW WWW 

Supersolid-2 AMI-2 

WWW WWW 

FIG. 4: (Color online) Schematic picture of the four different 
phases occurring in the Bose-Fermi mixture for bosonic filling 
3/2 and fermionic filling 1/2. We identify the Supersolid-1 
phase in which superfluidity coexists with a charge density 
wave with ANb < \. The AMI-I has localized bosons with 
AiVfc = |. The Supersolid-2 phase is defined by superfluidity 
coexisting with a charge density wave with i < ANb < |. 
The AMI-II has localized bosons with AN b = § . 



amplitude of the density oscillations as a function of the 
interspecies interaction Ufb- The amplitude of the den- 
sity oscillations is defined as ANf^ = \\n{ — n-^l, 
where ±1 refers to the two sublattices. The results show 
that the oscillation amplitude is a smooth function of 
Ufb for fermions and bosons. We identify four differ- 
ent regimes in the system. Schematic pictures for these 
four phases are given in Fig. 0J For weak interactions be- 
tween fermions and bosons the system is in the supersolid 
phase: the bosons are superfluid and there is a sponta- 
neous particle density oscillation in the system, which 
increases with increasing interaction Ufb- For some crit- 
ical Ufb the bosonic density amplitude reaches 1/2. At 
this point, the system undergoes a transition into the 
AMI-1 phase. Here the bosonic density is alternating 
between 1 and 2 on neighboring lattice sites. If we con- 
tinue to increase the interaction, only the amplitude of 
the fermionic density oscillations slowly increases. This 
continues up to the second phase transition from the AMI 
phase into second supersolid phase. In this region, with 
increasing Ufb, both amplitudes of the density oscilla- 
tions of fermions and bosons continuously increase, until 
the amplitude of the bosonic density oscillations reaches 
3/2. At this point a phase transition occurs from the 
supersolid into a second AMI phase. Within this AMI-2 
phase the bosons order themselves by alternating and 3 
bosons per site. Upon further increase in the interspecies 
interaction, the bosonic density oscillation - within our 
approximation - does not change, while the amplitude of 
the fermionic density oscillations converges to 1/2. In 
contrast to the case of half-filled hard-core bosons^, the 
superfluid order parameter is different on the two sub- 
lattices for this case, because there is no particle-hole 
symmetry for the bosons. This is visible in the inset of 
Fig. [3l where the superfluid order parameter on the two 
sublattices is plotted. 

An important observation concerns the order of the 
phase transitions. In the case of half-filled bosons, the 
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FIG. 5: (Color online) Spectral functions for the different 
phases. The parameters are chosen the same as in Fig. [3] 
The dashed green line corresponds to the supersolid- 1 phase 
(Ufb = GAD), the dash-dotted red line corresponds to the 
AMF1 phase with bosonic CDW oscillation 0.5 (Ufb = D) and 
the full blue line corresponds to the supersolid-2 phase (Ufb = 
1.95D). In the inset we plot the same spectral functions, at 
smaller frequencies. 

transition between the supersolid and AMI phase is a 
first order quantum phase transition^. However, for the 
bosonic density of 3/2 studied here, we find the transition 
to be of second order, as can be inferred from the lack of 
coexisting phases and the smooth behavior of all order 
parameters. 

We also study the local spectral functions in the differ- 
ent phases. The results are displayed in Fig. El The gap 
in the first supersolid phase is very small, as also found 
for the supersolid phase with half-filled bosons^. In the 
AMI phases we find that the fermions have a rather large 
gap. A more interesting structure emerges in the spectral 
function of the second supersolid phase. In this phase, in 
addition to the Hubbard sub-bands an additional peak 
arises in the spectral function. We have investigated the 
nature of the excitations responsible for this additional 
peak. These excitations correspond to a breaking of the 
alternating boson-fcrmion order in the system and there- 
fore indicate the instability of the system to phase sep- 
aration, which has only a slightly higher energy. In the 
AMI phase this energy difference is higher than in the 
supersolid phase, because the superfluid order parameter 
in the supersolid is oscillating (as seen from the inset of 
Fig. [3]) and therefore reduced. This leads to an increase 
in the energy and therefore enhances the instability to- 
wards phase separation. 

B. Strong coupling 

To gain a better analytic understanding of the sys- 
tem, we also consider a strong coupling approach. We 
propose a simple model, where in one of the sublattices 
on each site a fermion is localized, whereas the sites of 
the other sublattice are occupied by localized pairs of 



bosons. In addition we consider half-filled bosons on top 
of this arrangement. Within this model the AMI-l phase 
is described by the localization of the additional bosons 
on the "fermionic" sublattice. The AMI-2 phase corre- 
sponds to localization in the sublattice with the boson- 
pairs. The supersolid corresponds to the case where the 
additional bosons arc superfluid and dclocalized over all 
lattice sites. To describe the phase transition within this 
toy-model, we have to study localization of half-filled 
bosons in a superlattice. The effective Hamiltonian in 
the Gutzwillcr approach describing this situation has the 
form Ti-ef f = (W-i +Hi), where L is the number of 
lattices sites and 

Ht = -zt W -x (a\ + a x ) - (U b - ^p) (ni - \) (15) 

W_ 1= -zt b <p 1 V3(al 1 + a_i) + {U b - - ±) (16) 

where the index ±1 corresponds to the two sublattices. 
The sublattice marked by 1 is occupied by localized 
fermions and on each site of sublattice —1 there are two 
localized bosons. We have treated the additional boson as 
hard-core, which is justified because of the large bosonic 
interaction Ub- The factor -s/3 comes from the fact that 
in the second sublattice we have three bosons. We solve 
this system self-consistcntly and find the values when this 
system has a non-trivial solution (<p±i ^ 0). Our result 
shows that the system is superfluid in the following range: 

2U b - 2V3zt b < Ufb < 2U b + 2\fHzt b 

Also we compare the superfluid order parameter calcu- 
lated by strong coupling and GDMFT (see Fig. [5]). Our 
results show good agreement between these two results. 
Compared to the GDMFT- results, the strong coupling 
data are shifted towards smaller Bose-Fermi interaction. 
This shift is due to screening caused by the fact that in 
the superfluid phase the fermions are completely local- 
ized at the one sublattice, as wc assumed in this strong- 
coupling argument. In reality, due to virtual hopping 
processes, there is also a finite density of fermions on the 
other sublattice. This effectively reduces the interaction 
between fermions and bosons. 



V. SUMMARY 

We have investigated a Bose-Fermi mixture in an 
optical lattice by means of Generalized Dynamical 
Mean-Field Theory (GDMFT). This method consists of 
Gutzwiller mean-field for the bosons, and Dynamical 
Mean-Field Theory for the fermions, which are coupled 
onsite by the Bose-Fermi density-density interaction. We 
derived the self-consistency equations and showed that 
this method is well-controlled in the limit of high lattice 
coordination number z. 

We have applied the GDMFT scheme to a Bose-Fermi 
mixture with half-filled fermions, such that an instability 
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FIG. 6: (Color online) Superfluid order parameter on the two 
sublattices (a = ±1) as a function of the Fermi-Bose interac- 
tion, obtained by means of GDMFT and the strong coupling 
model. Parameters are chosen the same as in Fig. [3] In 
the inset we plot the same data, but the strong coupling re- 
sults are shifted towards stronger Ufb to compensate for the 
screening caused by virtual hopping processes of the fermions, 
which are not included in the toy-model. 

towards charge density-wave formation and hence super- 
solid order is present. We considered a bosonic filling of 
Nt, = 3/2, which allows for a series of phase transitions. 
A supersolid phase at small U fb is succeeded by an alter- 
nating Mott Insulator with alternating bosonic fillings 
1 and 2 for larger Ufb- For even larger Ufb & second 
supersolid phase is stable, untill for very large Ufb the 
ground state is formed by an AMI phase with alternating 
bosonic fillings and 3. The quantum phase transitions 
found here are of second order, in contrast to the case 
of half-filled bosons, where a first-order quantum phase 
transition was observed^. The phase diagram obtained 
here is particularly interesting because of the large am- 
plitude of the supersolid density oscillations between the 
two AMI phases, which will make experimental observa- 
tion easier. To compare quantitatively with experiments, 
it is necessary to perform calculations on the cubic lat- 
tice. This is beyond the scope of the current paper and 
will be pursued in the future. 
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APPENDIX A: DERIVATION OF THE 
EFFECTIVE ACTION 

To derive the self-consistency relations, we use the path 
integral formalism. The partition function of the Hamil- 
tonian ([1]) is given by : 

Z = JllD^Dc^D^Db.e- 3 (Al) 

The action is written as S = So + S° + AS, with 
S = [ dr\ V c* 0a {d T - fief ) c 0a + b* (d T - /x b ) b 

+Ufh^n{ [ + -jh b {n b - 1) + U fb h f h b Q } 
AS = - f dr{tfJ2' iStoh, + S&rSo*) (A2) 

S° = I dr{-t f J2 Kr^-tb^bfbj 

+ E (UfMA + Y»i( fi J - 1) + U fb n{n^ | 

where j3 is the inverse temperature, r is imaginary 
time, Ci,y, c* CT are the Grassmann variables describing the 
fermions, bi, b*, n\, hf are the usual C-numbers describ- 
ing the bosons and the number of fcrmions/bosons. Here 
the action is divided into three parts. So describes the 
"impurity site" , S° describes the system without the im- 
purity and AS is the coupling between them. means 
that the summations run only over the nearest neighbors 
of the "impurity site" and indicates a summation 

over all pairs of nearest neighbor sites excluding the "im- 
purity site" (i.e. i, j ^ 0). 

We now derive an effective action for the "impurity", 
defined by 




Using Eqs. EU [A2] and [A3] and with the definition 
AS = J drAS(T) we obtain 
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-So 



eff 



II D^Dc^DbtDb, e- s °e 



-AS 



S 



i^O.a 



Z 



J] D^Dc^DbtDl^Y, 



{-AS) n 



i#0,er 



n=0 



e ~ So ^[ 1 -J o dr(AS(T))° + - ^ dn J o dT 2 (AS(T 1 )AS(T 2 ))° + ...\ 

\ i i,j,cr 

-f 2 bl dT ^J Q bS(Ti)G^(ra - r 2 )b (r 2 ) + . 



(A4) 



where Z° is the statistical sum without the "impurity" site and (•■■)° are expectation values in the system not 
including the "impurity site". We have introduced the Nambu-spacc vector bo(r) = ( J*(r) J ' = ^ >i ^ T ^° as 

the bosonic supcrfluid parameter, G° ^(ti — t 2 ) = — (Tc^ti)^^)) as the Green's function for the fermions and 

G£ ^ (n - r 2 ) = -(r ^ ^ (SJfo^fo)) )° as the Green's function for the bosons in Nambu space. 

By the linked-cluster theorem we obtain 



Seff = So - h [ drY^mirKir) + c.c) + 1) £ / d n f dr 2 ^'^ (T (T 1 )G?. CT (n - r 2 )c 0rT (r 2 ) 

J a JO JO ida 

+hl fdr x f dT^Ki^Gl.^r, - r 2 )b (r 2 ) + . . . (A5) 
z Jo Jo 



I 

In this sum also higher order correlation functions ap- APPENDIX B: DERIVATION OF THE KINETIC 

pear (indicated by the dots). ENERGY 



In order to retain a finite kinetic energy, the hop- 
ping parameters should be rescalcd. The bosonic hop- 
ping parameter should be rescaled as if, = t^/z, and 
only the leading bosonic term describing the coupling 
to the bosonic superfluid order parameter survives in 
infinite dimensions. The fcrmionic hopping parameter 
will be rescalcd as tf = t^j^fz according to fermionic 
DMF T 49 ! 50 . After rescaling the hopping parameters and 
considering the limit z — > oo only the leading term for 
fermions and bosons survives. We obtain that Eq. IA5I 
reduces to the following relation: 



The fermionic kinetic energy is given by (to simplify 
the notations, we drop the summation over a): 



(Bl) 



where (ij) means summation over nearest neighbors. We 
now introduce the fermionic creation operators in the 
energy eigenbasis: 



(B2) 



S eff =S Q -t b [ dTj2'($°(T)b*(T) + C.c) (A6) 
Jo 

+*/£ / dT i / ^^ CT (n)G° > (T 1 -r 2 )c 0CT (T 2 ) 
Jo Jo 



where N is the number of lattice sites. The inverse trans- 
formation has the following form: 



(B3) 
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The following condition ensures that after the transfor- 
mation the Hamiltonian becomes diagonal: 



Using Eqs. ED E! E! EH iBTOl and IBTTI we obtain 



t 

'N 



EY V* = — (X X* ■+- X X* ) 

ni jn' at- / j \ ni jn' ' nj in' J 



-*E 



(>j) 



N/2 



+«>- w/2 / 1 \ 

= ~ E - ^«'MB4) = "* E E (^*A4A< -i + fc-cj 



+ < i i>- 



where a (ij)& denotes summation over the nearest neigh- 
bors such that i belongs to sublattice a and j belongs to 
sublatticc a = —a. At this point we have assumed that 
the lattice is bipartite. The second equality is based on 
the fact that both sublattices are identical and therefore 

For a bipartite lattice one can reverse the sign of the 
fermion creation/annihilation operators on one of the 
sublattices. This again yields an eigenstate of the Hamil- 
tonian (|B1[) . but with opposite sign. From this it directly 
follows that for each single-particle state with energy e n , 
there exists a state with energy — e„, i.e we can label the 
eigenstates such that 



~-n+N/2 — ~£n- 



(B5) 



From the Eqs. IB4I and |B51 it then follows that: 



N/2 

E 



2t 
N 



E X m X ]n' + h - C 



N/2 



(ij)- 

E £ " (4,1^,-1 + h.c 
n=l 

1 N 

o E £n (^n,l^n,-l + h 



(B12) 



In the last step we have used Eqs. IB5l and lB9l as follows: 

N/2 N/2 

E £ "^n,l^n,-l = E(~ £ "+A r /2)^j+jV/2.l(~^n+jV/2.-l) 



n=l 



n=l 



N 



El £ ™^n,l^n,-l 
n=JV/2+l 



XieSi,n+N/2 — Xi n and Xj e s^ 1 ,n+N/2 — —Xj n (B6) 



where S a (a = ±1) denotes the set of lattice points in 
sublattice a. 

Now we introduce two new operators 



£n,l ^ (Cn + C n+N/2 ) - -±= g X n 



(B7) 



Cn,-1 



1 

72 



Here and later we work modulo N, i.e. n + N = n. From 
Eqs. IB7l and [B8l onc easily obtains the following identity: 



C n +N/2,±1 — ±Cn,±l 



(B9) 



The inverse transformation has the following form: 



JV/2 

CieSi = ^jjp E X mC„,i (BIO) 

JV/2 

AjGS-i = E ( BU ) 



The next step is to go from summation to integral, 
and to take the expectation value of the kinetic energy 
operator. We obtain: 



Skin = \{j ds P°( £ ) £ (4,1^,-1 + h.cj^ 
1 f°° / 

= Urn - J de p°(e)e [(6^(0%^)) 
+ (cJ,_ 1 (0)c E , 1 (r))) 

/OO 
de p°{e)sB{e,T) 
- OO 

/OO 
de p°(e)e^T e- iuj " T B(e,uj n ) 
n 

/OO 
de p°(e)eB(e,Lu n ) 

n -°° 

/oo />oo 
dep°{e)e duf{uj)B{e,Lu + ), (B13) 
-oo J —oo 

where B( e ,r) = £ ((^(0)4,-1^) + <4,-i(0)4,i(t)>) 
and B = -±lmB 

TV 

These two terms are just the off-diagonal terms of the 
following Green's function matrix, which according to the 
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Dyson equation has the form 
G~ 1 (s,uj n ) ( lUJn + fif 



— E 
lLU n + /if 

Si(w) 


idln + /if — Si 

— £ 



(B14) 



We obtain 

G(e, w„) 

where 



C-i 



e C-i 



(B15) 



(B16) 



Therefore 
B(e,wn) 



CiC-i - £ 2 
1 / 1 



(B17) 



As one can easily see, the integral in Eq. IB13l stavs the 
same if we replace B(e,uj n ) by the following expression: 



1 



(B18) 



The advantage of this representation is that in the limit 
of one-sublattice it will reduce to the "usual" equation of 
the spectral function. 



APPENDIX C: DERIVATION OF THE 
SELF-ENERGY 

To derive the single-particle self-energy we use the 
equation of motion 



u>((A,B)) + (( 



H,A 



A, B 



(CI) 



and (. . .) denotes the usual expectation value. 

The Bose-Fermi Hamiltonian is given by Eq. [10l We 
use the following commutation relations: 



fa 



H, Cka 



- u fb fM - Uffjth -yy feCT c fe(7 (C3) 



k 



= -ekaCka - Vkafa ■ (C4) 

Here a = —a denotes the opposite spin state. 

First we will use the equation of motion for the case 
when A = f a and B = f a . Inserting the commutator 
relation (|C3[) in the equation of motion (|C1[) we obtain 

+»f)((Lfl)) -u fb {{fMJt)) - 

-Uf{{MLfl))-^V k a{{c k(T J a )) = 1 (C5) 



To calculate {(cka,fa)) we again use equation of motion 
(|C1|) . but in this case with A = Cka and B = f\. Using 
Eqs. El and [CD this yields 

(w - e k .) {{c ka JD) - Vk((L ft)) = • (C6) 
Equations. IC5l and I C6l then lead to 

(w + nf- A a (Lo)) G a {u>) -U fb F fba (Lo) -C/ / F/ / (w)=l 

(C7) 

where ((f = G a {u)) is the single-particle Green's 
function and A ff (w) = V£ / (u> — Ska) the hybridiza- 
tion function. We also define ({fflb ,/*)) = F/ 6 (w), 
((LflL Pa)) = Fj f H- Comparing Eq. E3to 

Gaiio)- 1 =W + Hf- KW-Vaifo). (C8) 

we finally obtain 



where r\ = + if A and B are both fermionic operators 
and rj = — otherwise. The notation ((. . .)) means 



((A,B)) 



i j dt e l 



A(t),B 



(C2) 



Ga(t0) 



Ga(u) 



(C9) 
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